##########################################################################
## Replication syntax for Malet & Hegewald (2025, JEPP):
## "The Changing Geography of Support for European Integration in the Shadow of the Ukraine War"
## Results manuscript mrp
##########################################################################

# R version: 4.5.1 
# Platform: x86_64-apple-darwin20

# here 1.0.1
# tidyverse 2.0.0
# autoMrP 1.0.6
# eurostat 4.0.0
# sf 1.0.21
# ggpubr 0.6.1

# 1. Load packages #######################################################
library(here)
library(tidyverse)
library(autoMrP)
library(eurostat)
library(sf)
library(ggpubr)

# 2. Basic setup #########################################################
rm(list=ls())
i_am("02_results_manuscript_mrp.R")

# 3. Import data ##########################################################
eurobarometer_census <- read.csv("eurobarometer.csv") 
eurobarometer_survey <- subset(eurobarometer_census, select = -PROP)

# 4. Figure 5 #############################################################
max_cores <- parallel::detectCores()
eurobarometer_survey$studyno1 <- as.character(eurobarometer_survey$studyno1)
eurobarometer_census$studyno1 <- as.character(eurobarometer_census$studyno1)

# 4.1 Defense #############################################################
eurobarometer_survey_def_pre <- eurobarometer_survey %>% 
  filter(studyno1=="7848") %>%
  select(CountryCode, GEO, pop_density, 
         Proposals_CommonDefence, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_census_def_pre <- eurobarometer_census %>% 
  filter(studyno1=="7848") %>%
  select(CountryCode, GEO, PROP, pop_density, 
         Proposals_CommonDefence, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_survey_def_post <- eurobarometer_survey %>% 
  filter(studyno1=="7888") %>%
  select(CountryCode, GEO, pop_density, 
         Proposals_CommonDefence, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_census_def_post <- eurobarometer_census %>% 
  filter(studyno1=="7888") %>%
  select(CountryCode, GEO, PROP, pop_density, 
         Proposals_CommonDefence, SEX, AGE, EDU) %>% 
  na.omit()

mrp_pre <- auto_MrP(
  y = "Proposals_CommonDefence",
  L1.x = c("SEX", "AGE", "EDU"),
  L2.x = c("pop_density"),
  L2.unit = "GEO",
  L2.reg = "CountryCode",
  bin.proportion = "PROP",
  survey = eurobarometer_survey_def_pre,
  census = eurobarometer_census_def_pre,
  ebma.size = 0,
  cores = max_cores,
  best.subset = FALSE,
  lasso = FALSE,
  pca = FALSE,
  gb = FALSE,
  svm = FALSE,
  mrp = TRUE
)

mrp_pre_df <- as.data.frame(mrp_pre$classifiers)
mrp_pre_df$prepost <- "pre"

mrp_post <- auto_MrP(
  y = "Proposals_CommonDefence",
  L1.x = c("SEX", "AGE", "EDU"),
  L2.x = c("pop_density"),
  L2.unit = "GEO",
  L2.reg = "CountryCode",
  bin.proportion = "PROP",
  survey = eurobarometer_survey_def_post,
  census = eurobarometer_census_def_post,
  ebma.size = 0,
  cores = max_cores,
  best.subset = FALSE,
  lasso = FALSE,
  pca = FALSE,
  gb = FALSE,
  svm = FALSE,
  mrp = TRUE
)

mrp_post_df <- as.data.frame(mrp_post$classifiers)
mrp_post_df$prepost <- "post"

mrp_df <- rbind(mrp_pre_df, mrp_post_df)

mrp_df_wide <- mrp_df %>%
  pivot_wider(
    names_from = prepost,
    values_from = mrp
  ) %>%
  mutate(change = (post - pre)*100)  

nuts1 <- get_eurostat_geospatial(nuts_level = 1, year=2021)
nuts2 <- get_eurostat_geospatial(nuts_level = 2, year=2021)
nuts3 <- get_eurostat_geospatial(nuts_level = 3, year=2021)
nuts_shapes <- bind_rows(nuts1, nuts2, nuts3)
nuts_shapes <- nuts_shapes %>% rename(GEO = NUTS_ID)

map <- mrp_df_wide %>% left_join(nuts_shapes, by = "GEO")
map <- st_as_sf(map)

defense <- ggplot(data = map) +
  geom_sf(aes(fill = change), color = "white", size = 0.2) +
  scale_fill_gradient2(
    low = "#4575b4",  
    mid = "#f7f7f7",  
    high = "#d73027", 
    midpoint = 0,    
    name = "\U0394 Support (in %)"
  ) +
  theme_minimal() + 
  labs(title = "Defence") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom"
  )

# 4.2 Foreign ############################################################
eurobarometer_survey_def_pre <- eurobarometer_survey %>% 
  filter(studyno1=="7848") %>%
  select(CountryCode, GEO, pop_density, 
         Proposals_CommonForeign, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_census_def_pre <- eurobarometer_census %>% 
  filter(studyno1=="7848") %>%
  select(CountryCode, GEO, PROP, pop_density, 
         Proposals_CommonForeign, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_survey_def_post <- eurobarometer_survey %>% 
  filter(studyno1=="7888") %>%
  select(CountryCode, GEO, pop_density, 
         Proposals_CommonForeign, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_census_def_post <- eurobarometer_census %>% 
  filter(studyno1=="7888") %>%
  select(CountryCode, GEO, PROP, pop_density, 
         Proposals_CommonForeign, SEX, AGE, EDU) %>% 
  na.omit()

mrp_pre <- auto_MrP(
  y = "Proposals_CommonForeign",
  L1.x = c("SEX", "AGE", "EDU"),
  L2.x = c("pop_density"),
  L2.unit = "GEO",
  L2.reg = "CountryCode",
  bin.proportion = "PROP",
  survey = eurobarometer_survey_def_pre,
  census = eurobarometer_census_def_pre,
  ebma.size = 0,
  cores = max_cores,
  best.subset = FALSE,
  lasso = FALSE,
  pca = FALSE,
  gb = FALSE,
  svm = FALSE,
  mrp = TRUE
)

mrp_pre_df <- as.data.frame(mrp_pre$classifiers)
mrp_pre_df$prepost <- "pre"

mrp_post <- auto_MrP(
  y = "Proposals_CommonForeign",
  L1.x = c("SEX", "AGE", "EDU"),
  L2.x = c("pop_density"),
  L2.unit = "GEO",
  L2.reg = "CountryCode",
  bin.proportion = "PROP",
  survey = eurobarometer_survey_def_post,
  census = eurobarometer_census_def_post,
  ebma.size = 0,
  cores = max_cores,
  best.subset = FALSE,
  lasso = FALSE,
  pca = FALSE,
  gb = FALSE,
  svm = FALSE,
  mrp = TRUE
)

mrp_post_df <- as.data.frame(mrp_post$classifiers)
mrp_post_df$prepost <- "post"

mrp_df <- rbind(mrp_pre_df, mrp_post_df)

mrp_df_wide <- mrp_df %>%
  pivot_wider(
    names_from = prepost,
    values_from = mrp
  ) %>%
  mutate(change = (post - pre)*100)  

nuts1 <- get_eurostat_geospatial(nuts_level = 1, year=2021)
nuts2 <- get_eurostat_geospatial(nuts_level = 2, year=2021)
nuts3 <- get_eurostat_geospatial(nuts_level = 3, year=2021)
nuts_shapes <- bind_rows(nuts1, nuts2, nuts3)
nuts_shapes <- nuts_shapes %>% rename(GEO = NUTS_ID)

map <- mrp_df_wide %>% left_join(nuts_shapes, by = "GEO")
map <- st_as_sf(map)

foreign <- ggplot(data = map) +
  geom_sf(aes(fill = change), color = "white", size = 0.2) +
  scale_fill_gradient2(
    low = "#4575b4",  
    mid = "#f7f7f7",  
    high = "#d73027", 
    midpoint = 0,    
    name = "\U0394 Support (in %)"
  ) +
  theme_minimal() + 
  labs(title = "Foreign") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom"
  )

# 4.3 Enlargement ########################################################
eurobarometer_survey_def_pre <- eurobarometer_survey %>% 
  filter(studyno1=="7848") %>%
  select(CountryCode, GEO, pop_density, 
         Proposals_EnlargeEU, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_census_def_pre <- eurobarometer_census %>% 
  filter(studyno1=="7848") %>%
  select(CountryCode, GEO, PROP, pop_density, 
         Proposals_EnlargeEU, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_survey_def_post <- eurobarometer_survey %>% 
  filter(studyno1=="7902") %>%
  select(CountryCode, GEO, pop_density, 
         Proposals_EnlargeEU, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_census_def_post <- eurobarometer_census %>% 
  filter(studyno1=="7902") %>%
  select(CountryCode, GEO, PROP, pop_density, 
         Proposals_EnlargeEU, SEX, AGE, EDU) %>% 
  na.omit()

mrp_pre <- auto_MrP(
  y = "Proposals_EnlargeEU",
  L1.x = c("SEX", "AGE", "EDU"),
  L2.x = c("pop_density"),
  L2.unit = "GEO",
  L2.reg = "CountryCode",
  bin.proportion = "PROP",
  survey = eurobarometer_survey_def_pre,
  census = eurobarometer_census_def_pre,
  ebma.size = 0,
  cores = max_cores,
  best.subset = FALSE,
  lasso = FALSE,
  pca = FALSE,
  gb = FALSE,
  svm = FALSE,
  mrp = TRUE
)

mrp_pre_df <- as.data.frame(mrp_pre$classifiers)
mrp_pre_df$prepost <- "pre"

mrp_post <- auto_MrP(
  y = "Proposals_EnlargeEU",
  L1.x = c("SEX", "AGE", "EDU"),
  L2.x = c("pop_density"),
  L2.unit = "GEO",
  L2.reg = "CountryCode",
  bin.proportion = "PROP",
  survey = eurobarometer_survey_def_post,
  census = eurobarometer_census_def_post,
  ebma.size = 0,
  cores = max_cores,
  best.subset = FALSE,
  lasso = FALSE,
  pca = FALSE,
  gb = FALSE,
  svm = FALSE,
  mrp = TRUE
)

mrp_post_df <- as.data.frame(mrp_post$classifiers)
mrp_post_df$prepost <- "post"

mrp_df <- rbind(mrp_pre_df, mrp_post_df)

mrp_df_wide <- mrp_df %>%
  pivot_wider(
    names_from = prepost,
    values_from = mrp
  ) %>%
  mutate(change = (post - pre)*100)  

nuts1 <- get_eurostat_geospatial(nuts_level = 1, year=2021)
nuts2 <- get_eurostat_geospatial(nuts_level = 2, year=2021)
nuts3 <- get_eurostat_geospatial(nuts_level = 3, year=2021)
nuts_shapes <- bind_rows(nuts1, nuts2, nuts3)
nuts_shapes <- nuts_shapes %>% rename(GEO = NUTS_ID)

map <- mrp_df_wide %>% left_join(nuts_shapes, by = "GEO")
map <- st_as_sf(map)

enlargement <- ggplot(data = map) +
  geom_sf(aes(fill = change), color = "white", size = 0.2) +
  scale_fill_gradient2(
    low = "#4575b4",  
    mid = "#f7f7f7",  
    high = "#d73027", 
    midpoint = 0,    
    name = "\U0394 Support (in %)"
  ) +
  theme_minimal() + 
  labs(title = "Enlargement") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom"
  )

# 4.4 Energy #############################################################
eurobarometer_survey_def_pre <- eurobarometer_survey %>% 
  filter(studyno1=="7848") %>%
  select(CountryCode, GEO, pop_density, 
         Proposals_CommonEnergy, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_census_def_pre <- eurobarometer_census %>% 
  filter(studyno1=="7848") %>%
  select(CountryCode, GEO, PROP, pop_density, 
         Proposals_CommonEnergy, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_survey_def_post <- eurobarometer_survey %>% 
  filter(studyno1=="7902") %>%
  select(CountryCode, GEO, pop_density, 
         Proposals_CommonEnergy, SEX, AGE, EDU) %>% 
  na.omit()

eurobarometer_census_def_post <- eurobarometer_census %>% 
  filter(studyno1=="7902") %>%
  select(CountryCode, GEO, PROP, pop_density, 
         Proposals_CommonEnergy, SEX, AGE, EDU) %>% 
  na.omit()

mrp_pre <- auto_MrP(
  y = "Proposals_CommonEnergy",
  L1.x = c("SEX", "AGE", "EDU"),
  L2.x = c("pop_density"),
  L2.unit = "GEO",
  L2.reg = "CountryCode",
  bin.proportion = "PROP",
  survey = eurobarometer_survey_def_pre,
  census = eurobarometer_census_def_pre,
  ebma.size = 0,
  cores = max_cores,
  best.subset = FALSE,
  lasso = FALSE,
  pca = FALSE,
  gb = FALSE,
  svm = FALSE,
  mrp = TRUE
)

mrp_pre_df <- as.data.frame(mrp_pre$classifiers)
mrp_pre_df$prepost <- "pre"

mrp_post <- auto_MrP(
  y = "Proposals_CommonEnergy",
  L1.x = c("SEX", "AGE", "EDU"),
  L2.x = c("pop_density"),
  L2.unit = "GEO",
  L2.reg = "CountryCode",
  bin.proportion = "PROP",
  survey = eurobarometer_survey_def_post,
  census = eurobarometer_census_def_post,
  ebma.size = 0,
  cores = max_cores,
  best.subset = FALSE,
  lasso = FALSE,
  pca = FALSE,
  gb = FALSE,
  svm = FALSE,
  mrp = TRUE
)

mrp_post_df <- as.data.frame(mrp_post$classifiers)
mrp_post_df$prepost <- "post"

mrp_df <- rbind(mrp_pre_df, mrp_post_df)

mrp_df_wide <- mrp_df %>%
  pivot_wider(
    names_from = prepost,
    values_from = mrp
  ) %>%
  mutate(change = (post - pre)*100)  

nuts1 <- get_eurostat_geospatial(nuts_level = 1, year=2021)
nuts2 <- get_eurostat_geospatial(nuts_level = 2, year=2021)
nuts3 <- get_eurostat_geospatial(nuts_level = 3, year=2021)
nuts_shapes <- bind_rows(nuts1, nuts2, nuts3)
nuts_shapes <- nuts_shapes %>% rename(GEO = NUTS_ID)

map <- mrp_df_wide %>% left_join(nuts_shapes, by = "GEO")
map <- st_as_sf(map)

energy <- ggplot(data = map) +
  geom_sf(aes(fill = change), color = "white", size = 0.2) +
  scale_fill_gradient2(
    low = "#4575b4",  
    mid = "#f7f7f7",  
    high = "#d73027", 
    midpoint = 0,    
    name = "\U0394 Support (in %)"
  ) +
  theme_minimal() + 
  labs(title = "Energy") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom"
  )

# 4.5 Combine maps #######################################################
Fig5 <- ggarrange(
  defense, foreign, enlargement, energy,
  ncol = 2, nrow = 2, 
  label.x = 0.02, label.y = 0.98, 
  align = "hv", 
  common.legend = TRUE, legend = "bottom" 
)

png(file = here("Figure5.png"), width = 8, height = 8, units = "in", res = 800)
Fig5
dev.off()